#############################################################
# Data management & analysis for
# How Do Voters Perceive Changes to the Rules of the Game? 
# Evidence from the 2014 Hungarian Elections
# Ahlquist, Ichino, Wittenberg & Ziblatt
# forthcoming, Journal of Comparative Economics
# Created on: 12/5/2017
# Created by: JSA
# Last updated on: 12/8/2017 
# Last updated by: NI 
#############################################################

library(MASS)
library(ggplot2)
library(foreign)
library(ri)
library(apsrtable)
library(multcomp)
library(RItools)
library(Matching)
library(lsmeans)
library(grid)
library(lattice)
library(stargazer)
library(MultinomialCI)
library(xtable)
rm(list=ls()); gc()

#helper functions
mcp.xtable<-function(x){
	require(xtable)
	pq<-x$test
	mtests <- cbind(pq$coefficients, pq$sigma, pq$tstat, pq$pvalues)
	error <- attr(pq$pvalues, "error")
	pname <- switch(x$alternativ, less = paste("Pr(<", ifelse(x$df ==0, "z", "t"), ")", sep = ""), 
		greater = paste("Pr(>", ifelse(x$df == 0, "z", "t"), ")", sep = ""), 
		two.sided = paste("Pr(>|",ifelse(x$df == 0, "z", "t"), "|)", sep = ""))
	colnames(mtests) <- c("Estimate", "Std. Error", ifelse(x$df ==0, "z value", "t value"), pname)
	xtable(mtests)}

vplayout <- function(x,y){viewport(layout.pos.row = x, layout.pos.col = y)}

gg_color_hue <- function(n) {
  hues = seq(15, 375, length=n+1)
  hcl(h=hues, l=65, c=100)[1:n]
} #emulate ggplot color hues

##settings
set.seed(456)  #note that ggplot automatically resets the seed in some graphics
setwd("YOUR DIRECTORY")
load("AIWZreplicationData.RData")
panel.w$fideszQuad <- relevel(panel.w$fideszQuad, ref="FideszDefector") #reference level
bothwaves<-subset(panel.w, panel.w$WAVE=="Both waves")  #subjects in both waves of the survey 

################################################################
## normalize (z transformation with the mean and sd of the controls only)
################################################################
m1 <- mean(panel.w$pre.legit[panel.w$treat=="cntrl"])
sd1 <- sd(panel.w$pre.legit[panel.w$treat=="cntrl"])
panel.w$pre.legit.norm <- (panel.w$pre.legit-m1)/sd1

m2 <- mean(na.omit(panel.w$post.legit[panel.w$treat=="cntrl"]))
sd2 <- sd(na.omit(panel.w$post.legit[panel.w$treat=="cntrl"]))
panel.w$post.legit.norm <- (panel.w$post.legit-m2)/sd2

m3 <- mean(na.omit(panel.w$d.legit[panel.w$treat=="cntrl"]))
sd3 <- sd(na.omit(panel.w$d.legit[panel.w$treat=="cntrl"]))
panel.w$d.legit.norm <- (panel.w$post.legit-m3)/sd3

################################################################
## summaries and values reported in the paper
################################################################

n.pre<-dim(panel.w)[1]
n.post<-sum(panel.w$WAVE=="Both waves")
## awareness of the electoral reforms prior to the election;
100*summary(panel.w$LSQ7)/n.pre  #89 % heard of reforms, pre-election

################################################################
## TABLE: Table 1
################################################################
panel.w$treat.num<-0
panel.w$treat.num[panel.w$treat=="Info"]<-1
panel.w$treat.num[panel.w$treat=="InfoParty"]<-2

## we are looking at $mean.Tr, $mean.Co, $tt$p.value for the 
## "pre-matching" dataset (since we did no matching) in the output below

## comparing control to "Info"
MatchBalance(treat.num~ fidesz + fideszQuad + ed + efficacy + age + gender + 
                        income + region + planned.turnout + LSQ23, 
                      data=panel.w[panel.w$treat.num!=2,])
## comparing control to "Info+Party"
MatchBalance(treat.num==2 ~ fidesz + fideszQuad + ed + efficacy + age + gender + 
                        income + region + planned.turnout + LSQ23, 
                      data=panel.w[panel.w$treat.num!=1,])
## comparing control to "Info" and "Info+Party" together
panel.w$InfoInfoParty <- panel.w$treat!="cntrl"
MatchBalance(InfoInfoParty ~ fidesz + fideszQuad + ed + efficacy + age + gender + 
                        income + region + planned.turnout + LSQ23, data=panel.w)

################################################################
### PLOT: Figure 1a
################################################################
pre.tab<-table(bothwaves$pre.legit,bothwaves$treat)
ci.pre<-apply(pre.tab,2,multinomialCI,alpha=.05)
ci.pre.cnt<-matrix(0,nrow=nrow(ci.pre), ncol=ncol(ci.pre))
for(i in 1:ncol(ci.pre)){
	ci.pre.cnt[,i]<-round(ci.pre[,i]*sum(pre.tab[,i]))
}
lb.pre<-c(ci.pre.cnt[1:5,1],ci.pre.cnt[1:5,2],ci.pre.cnt[1:5,3])
ub.pre<-c(ci.pre.cnt[6:10,1],ci.pre.cnt[6:10,2],ci.pre.cnt[6:10,3])
pre.dat<-data.frame(out = rep(-2:2,3),
	count=as.vector(pre.tab),
	treat=sort(rep(colnames(pre.tab),5)),
	ub = ub.pre,
	lb = lb.pre)

post.tab<-table(bothwaves$post.legit,bothwaves$treat)
ci.post<-apply(post.tab,2,multinomialCI,alpha=.05)
ci.post.cnt<-matrix(0,nrow=nrow(ci.post), ncol=ncol(ci.post))
for(i in 1:ncol(ci.post)){
	ci.post.cnt[,i]<-round(ci.post[,i]*sum(post.tab[,i]))
}
lb.post<-c(ci.post.cnt[1:5,1],ci.post.cnt[1:5,2],ci.post.cnt[1:5,3])
ub.post<-c(ci.post.cnt[6:10,1],ci.post.cnt[6:10,2],ci.post.cnt[6:10,3])
post.dat<-data.frame(out = rep(-2:2,3),
	count=as.vector(post.tab),
	treat=sort(rep(colnames(post.tab),5)),
	ub = ub.post,
	lb = lb.post)

#Fig 1a top
p1<-ggplot(pre.dat,  aes(x = out, y=count, fill= treat))
p1<-p1 + geom_bar(stat="identity", position="dodge", alpha=.7) + 
	labs(title = "Election reforms & legitimacy, pre-election") + 
	scale_x_continuous(name="Effect of reforms on legitimacy", breaks = c(-2, 0, 2), 
		labels = c("big effect for worse", "no effect", "big effect for better")) +
	scale_y_continuous(limits=c(0,275))+ 
	geom_errorbar(aes(ymin=lb, ymax=ub), position=position_dodge(.9), width=.3)+
	guides(fill = guide_legend(title = "treament")) + 
	theme(panel.background = element_blank(), legend.position="bottom")

#Fig 1a bottom
p2<-ggplot(post.dat,  aes(x = out, y=count, fill= treat))
p2<-p2 + geom_bar(stat="identity", position="dodge", alpha=.7) + 
	labs(title = "Post-election") + 
	scale_x_continuous(name="Effect of reforms on legitimacy", breaks = c(-2, 0, 2), 
		labels = c("big effect for worse", "no effect", "big effect for better")) +
	scale_y_continuous(limits=c(0,275))+ 
	geom_errorbar(aes(ymin=lb, ymax=ub), position=position_dodge(.9), width=.3)+
	#guides(fill = guide_legend(title = "treament")) + 
	theme(panel.background = element_blank(), legend.position="none")
pdf("Fig1a.pdf")
	grid.newpage()
	pushViewport(viewport(layout = grid.layout(2,1)))
	print(p1, vp = vplayout(1,1))
	print(p2, vp = vplayout(2,1))
dev.off()

################################################################
### PLOT: Figure 1b
################################################################
se.mean<-function(x,...){sd(x, na.rm=T)/sqrt(length(na.omit(x)))}
pl.m<-aggregate(panel.w$pre.legit, by=list(panel.w$treat), FUN=mean, na.rm=TRUE)
pl.m$period<-0
pl.se<-aggregate(panel.w$pre.legit, by=list(panel.w$treat), FUN=se.mean)

pol.m<-aggregate(panel.w$post.legit, by=list(panel.w$treat), FUN=mean, na.rm=TRUE)
pol.m$period<-1
pol.se<-aggregate(panel.w$post.legit, by=list(panel.w$treat), FUN=se.mean)
pol.m<-rbind(pol.m, pl.m)
pol.se<-rbind(pol.se, pl.se)

out<-data.frame(pol.m, pol.se[,2])
names(out)<-c("treatment", "legit","period","se")
out$upper<-out$legit+qnorm(0.975)*out$se
out$lower<-out$legit-qnorm(0.975)*out$se
p1<-ggplot(out ,aes(y=legit, x=period, color = treatment)) 
dodge <- position_dodge(width=0.5) 
pdf("figure1b.pdf")
p1 + geom_point(aes(shape=treatment),size=3, position=dodge) + #geom_line(aes(group=treatment), position=dodge) +
	geom_errorbar(aes(ymax = upper, ymin=lower), width=0.2, position=dodge) +
	labs(x=NULL, y="percieved electoral legitimacy" ) +
	scale_x_continuous(breaks = c( 0, 1), 
		labels = c("pre-election", "post-election")) +
	theme(panel.background = element_blank())
dev.off()
rm(out)

################################################################
## TABLE: Table 2
################################################################

table(panel.w$pre.legit,panel.w$treat)
### t-test
pl.aov<-aov(pre.legit~treat, data=panel.w)
aov.out <- glht(pl.aov, linfct = mcp(treat = "Tukey"), alternative = "less")
summary(aov.out)
diffmeans.legit.pre<-mcp.xtable(summary(aov.out, test=adjusted(type="holm")))


################################################################
## TABLE: Table 3
################################################################
#model 1
legit.pre.base <- lm(pre.legit ~ treat, data= panel.w) #eq 1 from PAP
#model 2
legit.pre.cov <- lm(pre.legit ~ treat + age + I(age^2) + gender + 
	income + region + planned.turnout, data=panel.w)   #eq 2 from PAP

legit.pre.base.norm <- update(legit.pre.base, pre.legit.norm ~.) #normalized version
## treatment effects in units of SD of control group
# treatInfo      -4.624e-02  4.551e-02  -1.016  0.30964   
# treatInfoParty -1.313e-01  4.547e-02  -2.887  0.00392 **
legit.pre.cov.norm <- update(legit.pre.cov, pre.legit.norm ~.) #normalized version
## treatment effects in units of SD of control group
# treatInfo             -5.732e-02  4.501e-02  -1.273 0.203004    
# treatInfoParty        -1.426e-01  4.497e-02  -3.172 0.001531 ** 

#model 3
legit.pre.fiquad.noint<-lm(pre.legit ~ treat + fideszQuad, data= panel.w)
#model 4
legit.pre.fiquad.cov.noint<-update(legit.pre.fiquad.noint, .~. + age + I(age^2) + gender +
	income + region + planned.turnout)
#model 5
legit.pre.fiquad<-lm(pre.legit ~ treat*fideszQuad, data= panel.w)
#model 6
legit.pre.fiquad.cov<-update(legit.pre.fiquad, .~. + age + I(age^2) + gender +
	income + region + planned.turnout)

######################################################################
## PLOT: Figure 2a
######################################################################

lpf.lsm<-lsmeans(legit.pre.fiquad, c("treat","fideszQuad"))
confint(as.glht(lpf.lsm),)
out<-confint(lpf.lsm, adjust="mvt") #multivariate t distribution for 6 comparisons
out$fideszQuad<-factor(out$fideszQuad, levels=c("NonFidesz","FideszDefector","FideszConvert", "StrongFidesz"))

p1<-ggplot(out ,aes(y=lsmean, x=fideszQuad, color = treat)) #this is figure for CPS
dodge <- position_dodge(width=0.5) 
pdf("fig2a.pdf")
p1 + geom_point(aes(shape=treat),size=3, position=dodge) +
	geom_errorbar(aes(ymax = upper.CL, ymin=lower.CL), width=0.2, position=dodge) +
	labs(x=NULL, y="effect of reforms on election legitimacy", 
		title = "Fidesz support & election legitimacy \npre-election" ) +
	scale_y_continuous(limits=c(-1.6,1.1))+
	scale_x_discrete(#breaks = c(0, 1,2,3), 
		labels = c("Non-supporters", "Defectors", "Converts", "Partisans")) +
	theme(panel.background = element_blank())
dev.off()

##########################################################################
## TABLE: multiple comparison of means
##########################################################################
legit.pre.fi2 <- lm(pre.legit ~ treatfidesz, data=panel.w)
lfp.out <- glht(legit.pre.fi2, 
	linfct = mcp(treatfidesz = c(
		"Info.0 - cntrl.0 =0", 
		"InfoParty.0 - cntrl.0 =0", 
		"InfoParty.0 - Info.0 =0", 
		"Info.1 - cntrl.1 =0", 
		"InfoParty.1 - cntrl.1 =0", 
		"InfoParty.1 - Info.1 =0")
	)
)
summary(lfp.out)


legit.pre.hte <- lm(pre.legit ~ treatfidesz + treateduc + treateff, data=panel.w)
lfp.out <- glht(legit.pre.fi2, 
	linfct = mcp(treatfidesz = c(
		"Info.0 - cntrl.0 =0", 
		"InfoParty.0 - cntrl.0 =0", 
		"InfoParty.0 - Info.0 =0", 
		"Info.1 - cntrl.1 =0", 
		"InfoParty.1 - cntrl.1 =0", 
		"InfoParty.1 - Info.1 =0")
	)
)
summary(lfp.out)

table(panel.w$fideszQuad, panel.w$treat)


##############################################################
## post-election
##############################################################

legit.prepost.base <- lm(d.legit ~ treat, data=panel.w)
legit.prepost.cov <- update(legit.prepost.base, . ~ . + age + I(age^2) +
	gender + income + region + turnout)
legit.prepost.fi <- lm(d.legit ~ treat*fidesz, data=panel.w)

######################################################################
## TABLE: Table 4
######################################################################

table(panel.w$post.legit,panel.w$treat)
pol.aov<-aov(post.legit~treat, data=panel.w)
aov.out <- glht(pol.aov, linfct = mcp(treat = "Tukey"), alternative = "less")
summary(aov.out)
diffmeans.legit.post<-mcp.xtable(summary(aov.out, test=adjusted(type="holm")))

#################################################################
## PLOT: Figure 3 
#################################################################

panel.w$fideszQuad2<-panel.w$fideszQuad
levels(panel.w$fideszQuad2)<-c("Nonsupporter", "Convert",  "Defector", "Partisan")
pdf("fig3.pdf")
dotplot(prop.table(table(panel.w$d.legit, panel.w$fideszQuad2,panel.w$treat),c(2,3)),
	main= "Within-subject opinion change",
	ylab = "proportion in attachment-treatment category",
	xlab = expression(paste(Delta, " election legitimacy")),
	par.settings = list(superpose.symbol = list(pch = c(16,2,0),
		col = c(gg_color_hue(3)[1], gg_color_hue(3)[2], gg_color_hue(3)[3]),
		cex=1.1),
	strip.background=list(col=grey(.9))), horizontal=F,
	auto.key=list(space="top", columns=3, title = "treatment", cex=.8),
	#par.strip.text=list(col="white", font=2)
 )
dev.off()


##############################################################
## TABLE: Table 5
##############################################################

#model 1
legit.post.base <- lm(post.legit ~ treat, data= panel.w)
#model 2
legit.post.cov <- lm(post.legit ~ treat + age + I(age^2) + gender + 
	income + region + planned.turnout, data=panel.w)

legit.post.base.norm <- update(legit.post.base, post.legit.norm ~.) #normalized version
legit.post.cov.norm <- update(legit.post.cov, post.legit.norm ~.) #normalized version
#model 3 
legit.post.fiquad.noint<-lm(post.legit ~ treat + as.factor(fideszQuad), data= panel.w)
#model 4 
legit.post.fiquad.cov.noint<-update(legit.post.fiquad.noint, . ~ . + age + I(age^2) + 
	gender + income + region + planned.turnout )
#model 5
legit.post.fiquad<-lm(post.legit ~ treat*as.factor(fideszQuad), data= panel.w)
#model 6
legit.post.fiquad.cov<-update(legit.post.fiquad, . ~ . + age + I(age^2) + 
	gender + income + region + planned.turnout )

##############################################################
## PLOT: Figure 2b
##############################################################
lpf.lsm<-lsmeans(legit.post.fiquad, c("treat","fideszQuad"))
confint(as.glht(lpf.lsm),)
out<-confint(lpf.lsm, adjust="mvt") #multivariate t distribution for 6 comparisons
out$fideszQuad<-factor(out$fideszQuad, levels=c("NonFidesz","FideszDefector","FideszConvert", "StrongFidesz"))

p1<-ggplot(out ,aes(y=lsmean, x=fideszQuad, color = treat)) 
dodge <- position_dodge(width=0.5) 
pdf("figure2b.pdf")
p1 + geom_point(aes(shape=treat),size=3, position=dodge) +
	geom_errorbar(aes(ymax = upper.CL, ymin=lower.CL), width=0.2, position=dodge) +
	labs(x=NULL, y="effect of reforms on election legitimacy", 
		title = "Fidesz support & election legitimacy \npost-election" ) +
	scale_y_continuous(limits=c(-1.6,1.1))+
	scale_x_discrete(#breaks = c(0, 1,2,3), 
		labels = c("Non-supporters", "Defectors", "Converts", "Partisans")) +
	theme(panel.background = element_blank())
dev.off()

##############################################################
##############################################################
## SUPPLEMENTAL MATERIALS 
##############################################################
##############################################################
## TABLE: Appendix D table 3
##############################################################

#model 1
legit.pre.ed.noint <- lm(pre.legit ~ treat + as.factor(ed), data= panel.w)
#model 2
legit.pre.ed.cov.noint <- update(legit.pre.ed.noint, .~. + age + I(age^2) +
	gender + income + region + planned.turnout)
#model 3
legit.pre.ed<-lm(pre.legit ~ treat*as.factor(ed), data= panel.w)
#model 4
legit.pre.ed.cov <- lm(pre.legit ~ treat*as.factor(ed) +
	age + I(age^2) + gender + income + region + planned.turnout, data=panel.w)
legit.pre.ed.norm <- update(legit.pre.ed, pre.legit.norm~.) #normalized version
legit.pre.ed.cov.norm <- update(legit.pre.ed.cov, pre.legit.norm~.) #normalized version
#model 5
legit.pre.ef.noint <- lm(pre.legit ~ treat +as.factor(efficacy), data= panel.w)
#model 6
legit.pre.ef.cov.noint <- update(legit.pre.ef.noint, .~. + age + I(age^2) +
	gender + income + region + planned.turnout)
#model 7
legit.pre.ef <- lm(pre.legit ~ treat*as.factor(efficacy), data= panel.w)
#model 8
legit.pre.ef.cov <- lm(pre.legit ~ treat*as.factor(efficacy) +
	age + I(age^2) + gender + income + region + planned.turnout, data= panel.w)

##############################################################
## TABLE: Appendix D table 4
##############################################################
#model 1
legit.post.ed.noint <- lm(post.legit ~ treat + as.factor(ed), data= panel.w)
#model 2
legit.post.ed.cov.noint <- update(legit.post.ed.noint, .~. + age + I(age^2) +
	gender + income + region + planned.turnout)
#model 3
legit.post.ed<-lm(post.legit ~ treat*as.factor(ed), data= panel.w)
#model 4
legit.post.ed.cov <- lm(post.legit ~ treat*as.factor(ed) +
	age + I(age^2) + gender + income + region + planned.turnout, data=panel.w)
legit.post.ed.norm <- update(legit.post.ed, post.legit.norm~.) #normalized version
legit.post.ed.cov.norm <- update(legit.post.ed.cov, post.legit.norm~.) #normalized version
#model 5
legit.post.ef.noint <- lm(post.legit ~ treat +as.factor(efficacy), data= panel.w)
#model 6
legit.post.ef.cov.noint <- update(legit.post.ef.noint, .~. + age + I(age^2) +
	gender + income + region + planned.turnout)
#model 7
legit.post.ef <- lm(post.legit ~ treat*as.factor(efficacy), data= panel.w)
#model 8
legit.post.ef.cov <- lm(post.legit ~ treat*as.factor(efficacy) +
	age + I(age^2) + gender + income + region + planned.turnout, data= panel.w)


##############################################################
## TABLE: Appendix E table 5
##############################################################

xtable(xBalance(panel.w$WAVE!="Only the first wave" ~ fidesz + ed + efficacy + age + 
                  gender + income + region + planned.turnout + LSQ23 + fidesz10 , 
                data=panel.w ),
                caption="Comparing pre-election wave-only ($n$ = 1500) and 
                both-wave ($n$ = 1500) standardized differences in means with $z$-scores" )


##############################################################
## TABLE: Appendix E table 6
##############################################################
#model 1
legit.pre.base.bwo<-update(legit.pre.base, data=bothwaves)
#model 2
legit.pre.cov.bwo<-update(legit.pre.cov, data=bothwaves)
#model 3
legit.pre.fiquad.noint.bwo<-update(legit.pre.fiquad.noint, data = bothwaves)
#model 4
legit.pre.fiquad.cov.noint.bwo<-update(legit.pre.fiquad.cov.noint, data = bothwaves)
#model 5
legit.pre.fiquad.bwo<-update(legit.pre.fiquad, data=bothwaves)
#model 6
legit.pre.fiquad.cov.bwo<-update(legit.pre.fiquad.cov, data=bothwaves)


##############################################################
## TABLE: Appendix E table 7
##############################################################
#model 1
legit.pre.ed.noint.bwo <- update(legit.pre.ed.noint, data= bothwaves)
#model 2
legit.pre.ed.cov.noint.bwo <- update(legit.pre.ed.cov.noint, data= bothwaves)
#model 3
legit.pre.ed.bwo<-update(legit.pre.ed, data= bothwaves)
#model 4
legit.pre.ed.cov.bwo <- update(legit.pre.ed.cov, data=bothwaves)
#model 5
legit.pre.ef.noint.bwo <- update(legit.pre.ef.noint, data= bothwaves)
#model 6
legit.pre.ef.cov.noint.bwo <- update(legit.pre.ef.cov.noint, data=bothwaves)
#model 7
legit.pre.ef.bwo <- update(legit.pre.ef, data= bothwaves)
#model 8
legit.pre.ef.cov.bwo <- update(legit.pre.ef.cov, data= bothwaves)


##############################################################
## TABLE: Appendix E table 8
##############################################################
#model OL1
legit.ol.base <- polr(ordered(pre.legit) ~ treat, data= panel.w, Hess=T)
#model OL2
legit.ol.cov <- polr(ordered(pre.legit) ~ treat + age + gender + 
	income + region + planned.turnout, data=panel.w, Hess=T)
#model OL3
legit.ol.fi <- polr(ordered(pre.legit) ~ treat*fideszQuad, data= panel.w, Hess=T)
#model OL4
legit.ol.fi.cov <- polr(ordered(pre.legit) ~ treat*fideszQuad +
	age + gender + income + region + planned.turnout, data= panel.w, Hess=T)

##############################################################
## TABLE: Appendix E table 9
##############################################################

#model OL1
legit.ol.ed <- polr(ordered(pre.legit) ~ treat*as.factor(ed), Hess=T, data= panel.w)
#model OL2
legit.ol.ed.cov <- polr(ordered(pre.legit) ~ treat*as.factor(ed) +
	age + gender + income + region + planned.turnout, data=panel.w, Hess=T)
#model OL3
legit.ol.ef <- polr(ordered(pre.legit) ~ treat*as.factor(efficacy), data= panel.w, Hess=T)
#model OL4
legit.ol.ef.cov <- polr(ordered(pre.legit) ~ treat*as.factor(efficacy) +
	age + gender + income + region + planned.turnout, data= panel.w, Hess=T)


##############################################################
## TABLE: Appendix E table 10
##############################################################
## using only those who "heard of" reforms
informed<-subset(panel.w, panel.w$LSQ7=="Yes ")
uninformed<-subset(panel.w, panel.w$LSQ7!="Yes ")

#model 1
legit.pre.base.inf<-update(legit.pre.base, data=informed)
#model 2
legit.pre.cov.inf<-update(legit.pre.cov, data=informed)
#model 3
legit.pre.fiquad.noint.inf<-update(legit.pre.fiquad.noint, data = informed)
#model 4
legit.pre.fiquad.cov.noint.inf<-update(legit.pre.fiquad.cov.noint, data = informed)
#model 5
legit.pre.fiquad.inf<-update(legit.pre.fiquad, data=informed)
#model 6
legit.pre.fiquad.cov.inf<-update(legit.pre.fiquad.cov, data=informed)


##############################################################
## TABLE: Appendix E table 11
##############################################################

#model 1
legit.pre.ed.noint.inf <- update(legit.pre.ed.noint, data= informed)
#model 2
legit.pre.ed.cov.noint.inf <- update(legit.pre.ed.cov.noint, data= informed)
#model 3
legit.pre.ed.inf<-update(legit.pre.ed, data= informed)
#model 4
legit.pre.ed.cov.inf <- update(legit.pre.ed.cov, data=informed)
#model 5
legit.pre.ef.noint.inf <- update(legit.pre.ef.noint, data= informed)
#model 6
legit.pre.ef.cov.noint.inf <- update(legit.pre.ef.cov.noint, data=informed)
#model 7
legit.pre.ef.inf <- update(legit.pre.ef, data= informed)
#model 8
legit.pre.ef.cov.inf <- update(legit.pre.ef.cov, data= informed)


##############################################################
## TABLE: Appendix F table 12
##############################################################

set.seed(20150207)
panel.w$treatNum<-as.numeric(panel.w$treat)-1

## generate null distribution for difference-in-means statistic, without moderators
AIWZri<-function(Y,Z,level1,level0,numiter){
  dY.perm<-rep(NA,numiter)
  for(i in 1:numiter){
    s<-sample(Z)
    dY.perm[i]<-mean(Y[s==level1])-mean(Y[s==level0])
  }
  dY.perm
}

## generate null distribution for difference in difference-in-means statistic across levels of moderator X (heterogeneous treatment effects)
AIWZriX<-function(Y,Z,level1,level0,numiter,X,Xlevel1,Xlevel0){
  dY.perm<-rep(NA,numiter)
  for(i in 1:numiter){
    s<-sample(Z)
    dY.perm[i]<-mean(Y[s==level1&X==Xlevel1])-mean(Y[s==level0&X==Xlevel1])-(mean(Y[s==level1&X==Xlevel0])-mean(Y[s==level0&X==Xlevel0]))
  }
  dY.perm
}

Z<-panel.w$treatNum
numiter<-10000

Y<-panel.w$pre.legit

# hypothesis 1a
d1a<-mean(Y[Z==1])-mean(Y[Z==0])
r1a<-AIWZri(Y, Z, 1,0, numiter)
p1a<-mean(r1a<d1a)
p1a

# hypothesis 1b
d1b<-mean(Y[Z==2])-mean(Y[Z==0])
r1b<-AIWZri(Y, Z, 2,0, numiter)
p1b<-mean(r1b<d1b)
p1b

# hypothesis 1c -- two-sided test
d1c<-mean(Y[Z==2])-mean(Y[Z==1])
r1c<-AIWZri(Y, Z, 2,1, numiter)
p1c<-mean(abs(r1c)>abs(d1c))
p1c


# for family 4:

Y<-panel.w$pre.dem
Z<-panel.w$treatNum

# hypothesis 4a
d4a<-mean(Y[Z==1])-mean(Y[Z==0])
r4a<-AIWZri(Y, Z, 1,0, numiter)
p4a<-mean(r4a<d4a)
p4a

# hypothesis 4b
d4b<-mean(Y[Z==2])-mean(Y[Z==0])
r4b<-AIWZri(Y, Z, 2,0, numiter)
p4b<-mean(r4b<d4b)
p4b

# hypothesis 4c -- two-sided test
d4c<-mean(Y[Z==2])-mean(Y[Z==1])
r4c<-AIWZri(Y, Z, 2,1, numiter)
p4c<-mean(abs(r4c)>abs(d4c))
p4c

# hypothesis 6b
X<-as.numeric(panel.w$fidesz)-1
d6b<-mean(Y[Z==2&X==0])-mean(Y[Z==1&X==0])
r6b<-rep(NA,numiter)
for(i in 1:numiter){
  s<-sample(Z)
  r6b[i]<-mean(Y[s==2&X==0])-mean(Y[s==1&X==0])
}
p6b<-mean(r6b<d6b)
p6b


# hypothesis 7a
Y<-panel.w$d.legit
df<-na.omit(data.frame(Y,Z))
Y<-df$Y
Z<-df$Z

d7a<-mean(Y[Z==1])-mean(Y[Z==0])
r7a<-AIWZri(Y,Z,1,0,numiter)
p7a<-mean(r7a<d7a)
p7a

# hypothesis 7b
d7b<-mean(Y[Z==2])-mean(Y[Z==0])
r7b<-AIWZri(Y,Z,2,0,numiter)
p7b<-mean(r7b<d7b)
p7b

# hypothesis 7c -- two sided
d7c<-mean(Y[Z==2])-mean(Y[Z==1])
r7c<-AIWZri(Y,Z,2,1,numiter)
p7c<-mean(abs(r7c)<abs(d7c))
p7c

###END